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Introduction 

This  is  a  report  on  the  work  done  during  the  third  year  of  of  this  grant  plus  a  six 
month  no-cost  extension.  The  project  was  to  develop  theoretical  models  of  several  aspects 
of  accelerating  arc  plasmas  in  electromagnetic  railguns  with  plasma  driven  projectiles. 

The  grant  provided  three  months  of  support  during  the  summer  of  1986  for  Dr.  Huerta, 
and  a  further  two  months  during  the  summer  of  1987.  Dr.  J.  C.  Nearing  was  supported  for 
two  months  during  the  summer  of  1986.  The  other  person  who  received  support  was  Mr. 
G.  Christopher  Boynton  who  is  a  graduate  student  doing  his  dissertation  on  this  work.  He 
was  supported  fully  during  the  whole  year  and  a  half. 

A  group  of  researchers  on  railgun  problems  was  developed  at  the  University  of  Miami 
largely  due  to  the  support  flowing  from  this  project.  There  appeared  a  strong  possibility 
of  funding  by  the  SDIO(IST)  office  of  a  consortium  proposed  by  Dr.  Huerta,  Dr.  Harry  S. 
Robertson,  and  other  faculty  at  the  University  of  Miami.  Several  of  us  attended  the  23  June 
1986  EML  Science  Interchange  Workshop  at  the  Air  Force  Armament  Laboratory  of  Eglin 
Air  Force  Base  that  was  sponsored  by  Dr.  Leonard  Caveny  of  SDIO(IST).  There  I  presented 
our  research  plans  on  fundamental  processes  in  rapidly  moving  high  current  arcs.  Soon 
thereafter  we  received  two  letters  of  intent  to  fund.  Delays  kept  coming  up,  however.  Col 
David  Finkleman  asked  us  several  times  to  submit  extensively  revised  proposals  following 
SDIO/IST  instructions.  I  made  a  presentation  to  Drs.  Dwight  Dustin  and  Caveny,  at 
SDIO(ISDT)  headquarters  in  Washington  in  February,  1987.  Finally  our  proposal  was 
declined  declined  in  April  of  1987.  A  great  deal  of  time  and  energy  was  wasted  by  the  P. 
I.  and  others  during  this  process.  The  whole  episode  was  a  fiasco. 

Description  of  the  Technical  Work 

During  the  second  year  of  this  project  Miss  Ann  Decker  obtained  her  Ph.  D.  submitting 
a  dissertation  on  the  effects  of  finite  conductivity  and  compressibility  on  the  development  of 
the  Kruskal-Schwarzchild  instability.  The  inclusion  of  finite  conductivity  effects  produced 
a  fourth  order  eigenvalue  equation  for  the  modes.  She  did  a  numerical  solution  of  the 
equation  and  obtained  a  dispersion  relation.  Part  of  that  work  was  presented  by  the  P.  I. 
at  the  First  EM  Gun  Armature  Workshop,  24-26  June,  1986.  In  the  course  of  preparing 
this  work  for  publication  the  P.  I.  found  that  Miss  Decker  had  neglected  to  examine  the 
effects  of  a  singularity  in  the  fourth  derivative  term  of  the  mode  equation.  This  has  called 
for  a  revision  of  the  problem.  A  physicist  on  a  sabbatical  visit  from  Argentina,  Dr.  Felix 
Rodriguez-Trelles,  has  been  of  great  help  in  this. 

The  work  with  Dr.  Nearing  concentrated  on  modeling  the  current  in  the  rails  for 
arbitrarily  time  dependent  current  and  velocity.  Also  the  temperature  distribution  was 
obtained  in  important  cases.  A  paper  on  this  work  was  submitted  for  the  4th  Symposium 
on  Electromagnetic  Launch  Technology  in  Austin  Texas  during  April  11-14  1988,  and  for 
publication  in  the  IEEE  Transactions  on  Magnetics.  A  copy  of  this  paper  is  attached. 

The  work  with  Mr.  Boynton  centers  on  a  two  dimensional  time  dependent  simulation 
of  plasma  armatures.  We  use  the  equations  of  resistive  MHD  and  we  use  a  two  dimensional 
FCT  code  to  advance  all  quantities  in  time.  A  good  deal  of  effort  has  also  expended  on 
developing  graphical  methods  that  allow  good  display  of  the  results.  We  have  done  several 
runs  in  the  University  of  Miami’s  VAX  8650,  which  does  about  1,000  steps  per  hour  of 


CPU  time.  A  paper  with  some  of  the  results  of  this  work  for  an  adiabatic  case  was  also 
submitted  for  the  4th  Symposium  on  Electromagnetic  Launch  Technology  in  Austin  Texas, 
and  for  publication  in  the  IEEE  Transactions  on  Magnetics.  A  copy  of  this  paper  also  is 
attached.  We  have  results  beyond  those  prepared  for  the  Austin  conference.  The  last 
figures  attached  represent  the  pressure  in  the  plasma  developing  in  time  from  120,000  time 
steps,  87/isec  since  firing,  with  the  armature  4.6  cm  down  the  barrel,  in  five  figures  until 
the  last  one  at  89.87/xsec  with  a  travel  of  4.9  cm.  The  plasma  was  initially  in  mechanical 
equilibrium  with  a  2%  perturbation  of  the  pressure  near  the  rear  of  the  armature.  By  the 
time  the  armature  has  advanced  only  about  5  cm  serious  disruptions  appear  in  it.  Due  to 
the  violent  plasma  motions  the  time  step  has  to  be  reduced  as  we  go  along  to  satisfy  the 
Courant  condition  and  prevent  crashes  due  to  numerical  instabilities.  This  is  having  a  bad 
effect  on  the  CPU  time  requirement. 

Clearly  this  problem  is  one  that  requires  a  supercomputer.  During  the  summer  of 
1988  Dr.  Huerta  and  Mr.  Boynton  are  going  to  be  working  at  Eglin  AFB  under  the  USAF 
Summer  Faculty  Research  Program.  We  expect  to  be  able  to  use  the  Cyber  205  at  AFATL. 
This  would  be  of  great  benefit  to  our  research. 
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SKIN  AND  HEATING  EFFECTS  OF  RAILGUN  CURRENT 
J.  C.  Nearing  and  M.  A.  Huerta 
University  of  Miami 
Physics  Department 
Coral  Gables,  FL  33124 

Abstract 

Wc  present  a  calculation  in  a  simplified  geometry  for  the  current  distribution  in  the 
rails,  taking  into  account  the  motion  of  the  armature  and  the  time  variation  of  the  current. 
Closed  form,  asymptotic,  results  for  the  current  density  are  obtained  for  arbitrary  time 
dependent  currents  and  velocities,  in  the  limit  in  which  the  length  scale  ik^1  is  small,  where 
ko  =  a  is  the  electrical  conductivity  of  the  rails,  and  v  is  the  speed  of  the  armature. 

Because  of  eddy  current  effects  the  rail  current  may  reverse  in  portions  of  the  rails  when 
the  total  current  decreases.  The  current  is  used  as  a  source  of  Joule  heating  to  find  the 
temperature  distribution  in  the  rails.  The  heat  diffusivity  is  negligible  and  we  are  able  to 
give  numerical  results  concerning  melting.* 

Introduction 

In  a  railgun  the  current  goes  into  one  rail,  passes  through  the  armature,  and  returns  via 
the  other  rail.  The  armature  is  immersed  in  the  magnetic  field  produced  by  the  rail  currents 
and  experiences  a  strong  magnetic  force  that  pushes  the  projectile  in  front  of  it.  The 
distribution  of  current  in  the  rails  is  one  of  the  important  factors  in  the  operation  of  railguns 
because  it  determines  the  rail  resistance  and  the  joule  heating  losses.  A  great  deal  of  work 
has  been  done  to  model  this  current  distribution.  KerriskJ’2  did  a  numerical  calculation 
of  the  current  distribution  and  of  the  accompanying  temperature  distribution,  allowing 
the  electrical  conductivity  to  depend  on  temperature  and  the  magnetic  permeability  to 
depend  on  magnetic  field  strength.  He  treated  carefully  the  two  dimensional  variations 
in  the  rectangular  cross  section  of  the  rails  but  neglected  the  variations  in  the  direction 
along  the  rails.  Marshall3  discussed  qualitatively  the  problem  of  the  variation  along  the 
rails.  Long4  attempted  a  solution  for  a  case  with  steady  current  and  speed.  Drake  and 
Rathmann5  obtained  infinite  series  solutions  that  described  the  variation  of  the  skin  depth 
along  the  rail  due  to  the  motion  of  the  armature. 

In  a  railgun,  the  part  of  the  rails  ahead  of  the  armature  is  in  a  region  where  there  is 
almost  no  magnetic  field.  As  the  armature  sweeps  along  the  rails,  the  rails  are  exposed 
to  the  strong  field  behind  the  armature.  Due  to  the  rapid  motion  of  the  armature 
along  the  rails,  the  current  and  magnetic  field  do  not  have  time  to  diffuse  completely 
into  the  rails  and  are  concentrated  in  a  skin  layer  near  the  rail  surfaces.  We  seek  to 
describe  analytically  how  the  magnetic  field  and  current  diffuse  into  the  rails  and  what 
the  temperature  distribution  caused  by  the  Joule  heating  is.  The  important  skin  depth 
that  arises  can  be  written  as  8  =  where  t  is  the  time  since  the  magnetic  field 

started  to  diffuse  into  the  conductor  of  conductivity  a,  and  magnetic  permeability  hq. 
There  are  several  related  problems  presented  by  Knoepfel6  with  known  analytic  solutions. 

*  This  work  supported  in  part  by  the  Air  Force  Office  of  Scientific  Research  under  grant 
number  84-01 16. 


They  involve  the  calculation  of  the  magnetic  field  and  current  inside  a  semi-infinite  space 
bounded  by  a  plane,  or  inside  a  cylinder,  when  the  magnetic  field  outside  is  uniform  in 
space  and  has  a  known  time  dependence. 

In  order  to  elucidate  the  effects  of  the  motion  of  the  armature,  our  analytic  solution 
takes  into  account  the  variations  along  the  rail  direction,  x  as  shown  in  Fig.  1.  The  rails 
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are  taken  as  infinitely  thick  because  their  thickness  is  much  greater  than  the  skin  depth. 
We  also  take  into  account  variations  in  the  current  density  in  the  direction  into  the  rails, 
y.  This  will  enable  us  to  describe  skin  depth  effects  on  the  inside  surface  of  the  rails. 
We  simplify  the  geometry  by  neglecting  variations  in  the  z  direction,  in  effect  making  the 
problem  that  of  two  rails  of  infinite  height.  In  a  complete  treatment  of  the  problem  Jy(t) 
in  the  armature  would  be  determined  by  the  properties  of  the  armature  in  conjunction 
with  the  properties  of  the  rails.  We  model  the  armature  by  an  assumed  time  dependent 
distribution  of  current  density  in  an  armature  moving  to  the  right  along  the  rails  with  an 
arbitrary  velocity  v0(t)  =  dx0/dt. 

We  obtain  rigorous  general  results  for  the  current  distribution  in  the  rails  and  find  the 
temperature  distribution  for  special  cases.  We  get  conditions  for  melting  depending  on  the 
current,  and  much  less  strongly  on  the  length  of  the  arc.  We  also  give  an  expression  for 
the  breech  voltage  in  terms  of  the  current,  and  discuss  the  circumstances  when  it  reduces 
to  the  simple  expressions  that  are  typically  used  in  railgun  circuit  models.  One  of  the 
most  interesting  results  is  the  reversal  of  current  density  in  the  rails,  and  possibly  in  the 
armature,  whenever  the  total  rail  current  decreases. 

Current  Calculation 

In  the  coordinate  system  shown  in  Fig.  1,  the  upper  rail  occupies  the  space  y  >  0. 
The  numbers  1-6  are  used  later  to  identify  important  voltages.  The  lower  rail  occupies 
the  space  y  <  —w  and  carries  a  current  that  is  a  mirror  image  of  the  current  in  the  upper 
rail.  The  armature  between  the  rails  occupies  the  space  —w  <  y  <  0  and  has  a  current 
density  Jy(x,t) 

Jy(x,  t)  =  -I'(t)f(x  -  x0(t)) 


where  /'  is  the  current  per  rail  height,  /  represents  the  spatial  distribution  of  the  current 
in  the  armature,  and  x0(t)  is  the  position  of  the  center  of  the  armature.  /  is  normalized 


to  /  f(x)dx  =  1.  Everything  is  independent  of  z,  and  only  the  x,  y,  and  t  dependencies 
are  to  be  determined.  To  analyze  the  fields  resulting  from  this  source,  we  write  Jy  as  a 
Fourier  integral 


•MM)  =  J  /ffs( .  (-«.  <  y  <  0)  (1) 

/(*)  =  /  5?  9<*)ei‘ '•  M*.")  =  -/ 

In  Maxwell’s  equations  we  neglect  the  displacement  current  compared  to  the  conduc¬ 
tion  current  because  all  the  time  constants  considered  are  much  larger  than  e0/cr.  The 
resulting  equations  are 


__  „  .dB 

VxE  =  "5T 


V  x  B  =  AioJ  and 
In  the  domain  y  <  0,  a  vacuum,  the  magnetic  field  due  to  a  current  density  Jy  = 


is 


^  e*(fc2!-a,t)  (y  <  0). 

K 

Taking  the  full  Jj,  of  Eq.  (1)  into  account  we  find  the  magnetic  field  for  y  <  0  to  be 

/OO 

dx1  f{x')  (y  <  0),  (2) 

where  £  =  x  —  xo(t).  The  absolute  value  |£|  is  the  distance  along  the  rail  from  the  point 
with  coordinate  x  to  the  center  of  the  armature.  In  arriving  at  (2)  we  have  used  the 
integral 

f+°°  dkeikz  =  f  1  x<0 

2tt  7-00  k  1  0  x  >  0 


where  the  path  of  integration  goes  above  the  pole  at  k  =  0.  The  magnetic  field  of  Eq.  (2) 
is  constant  behind  the  armature,  where  £  <  0,  and  zero  ahead  of  it,  where  £  >  0,  as  one 
would  expect  in  this  geometry. 

In  the  domain  y  >  0,  the  interior  of  the  upper  rail,  the  electric  field  is  eliminated  using 
the  conductivity  equation,  J  =  <rE.  Eliminating  E  and  J  from  Maxwell’s  equations  gives 
the  diffusion  equation  for  the  magnetic  field 


V2B  =  n$(T 


3B 

at 


(3) 


The  vector  B  has  only  a  z-component.  A  single  Fourier  component  is  written  Bz  — 
F(y)  exp  [z(fcx  —  w<)]  (y  >  0).  When  this  is  substituted  in  Eq.  (3),  the  resulting  equation 


is 


dy  F  -f  (io7<r/io  -  k2)F  =  0, 
3 


whose  solution  is 


F{y)  =  eay ,  .  =  (4) 

The  coefficient  for  F(y )  comes  from  the  continuity  of  Bz  at  y  =  0. 

Since  the  fields  and  currents  must  die  off  as  y  becomes  large  and  positive,  the  sign  of 
the  square  root  is  3?(a)  <  0.  In  order  to  obtain  the  value  of  the  magnetic  field  for  y  >  0, 
we  integrate  over  k  and  u>  in  Eq.  (1).  The  value  of  B  for  y  >  0  will  then  be 

_  ,  .  .  [  ,  .  f  dk  [  du>  1 

l)  =  -tfi0J  dt  J  -J  -- 

The  a?  integral  can  be  done  with  no  further  assumptions  on  the  time  dependence  of 
or  the  shape  /  of  the  armature’s  current  distribution.  The  integral  J  eay~tw(t~t  '>du>/2‘n  is 
zero  for  t  <  t1 .  The  integrand  has  a  branch  point  at  w  =  -•-ik2  /  yo<r,  so  when  t  >  t\  we 
deform  the  contour  around  this  branch  toward  —zoo.  Change  variables  to  s  —  iw  —  fc2//i0o' 
and  observe  that  the  imaginary  part  of  the  integral  vanishes.  Using  integration  by  parts 
we  obtain 


(5) 


Bz  =  -ifi o  J '  dt'  /'(<')  J  ^  (6) 

V  /.  ~  -fea(t-t')/Mo<T  -yaMQg/4(t-t')  >  q 

2  y  7r(<  —  t')3  6  ,y>Ul 

The  integral  of  Eq.  (6)  is  a  general  expression  for  the  magnetic  field  in  this  geometry.  (The 
behavior  of  this  expression  as  y  — >  0  needs  to  be  treated  carefully.)  In  order  to  make  further 
progress  on  this  integral,  it  is  necessary  to  assume  a  shape  for  the  current  distribution  in 
the  armature.  We  have  considered  several  forms  for  f(x).  One  of  the  simplest  forms  to 
take  that  yields  interesting  results  is 


{-L<x<  L). 


With  this  form,  the  k  integral  of  Eq.  (6)  can  be  done.  It  is  of  the  form 


f  dk  1  sin  kL ^ikx_ak2 

J  2tt  k  kL 


The  contour  goes  from  — oo  to  +oo,  passing  above  the  singularity  at  k  =  0.  Bring  the 
contour  down  to  the  real-fc  axis,  and  we  get  a  principal  value  integral  plus  a  semi-circular 
contour  just  over  the  pole  at  zero.  Only  the  cosine  part  of  exp (ikx)  contributes  to  the 
latter,  and  only  the  sine  part  to  the  former.  This  integral  is  then 


We  combine  the  two  sines  into  the  difference  of  cosines  and  use  a  tabulated7  integral  to 
obtain 


Bz(x,y,t)  =  -noy]/^  fdt'nm-tr^ 

j^«  -  *0  -  A)  erf  ^(x  -  x0  -  L)y/n0o-/4(t  -  t') 
—  (x  —  x0  -f  A)  erf  J^(x  -  x0  +  L)y//j.o<r/4(t  —  V) 
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The  above  expression  for  the  magnetic  field  allows  an  arbitrary  thickness  L  for  the 
armature,  an  arbitraxy  time  dependence  for  the  current  I(t ),  and  a  general  time  dependence 
for  the  armature  position  xo(t).  The  curl  of  this  expression  for  Bz  will  give  the  current 
density.  Rather  than  writing  down  the  complicated  expressions  that  result  in  the  general 
case,  we  will  specialize  to  two  cases  of  interest.  In  the  first  case  we  consider  a  thin  armature 
( L  — ►  0)  with  a  general  time  dependent  current  and  velocity.  The  thin  armature  case 
simplifies  the  derivation  of  an  approximate  current -volt age  relation  for  the  railgun.  It 
also  allows  a  closed  form  analytical  expression  for  the  temperature  rise  caused  by  very 
concentrated  currents.  In  the  second  case  we  will  consider  a  thick  armature  but  with  the 
current  and  the  velocity  kept  constant  in  time.  This  case  is  of  interest  because  we  are  able 
to  calculate  the  temperature  rise  of  the  rails  for  realistically  large  arcs. 

Time  dependent  results  for  a  thin  armature 

A  thin  armature  could  represent  the  sort  of  current  spots  ,or  filaments,  that  are 
observed.  In  this  section,  then,  I'(i)  could  represent  the  current  in  a  spot,  and  not  the 
total  current  through  the  rail.  Of  course,  due  to  the  two  dimensional  nature  of  our  model, 
the  spot  is  really  modeled  as  a  sheet.  The  thin  armature  result  is  obtained  by  taking  the 
limit  as  L  — >  0  in  Eq.  (7) 

Bz{x,y,t)  =  J  dt' r(t'){t~t')~3/2  (8) 

.  e-y  fiocr/Mt-t  )  j’j  _  erf  s/]x0(rf4{t  -  U)  . 


The  curl  of  this  expression  for  Bz  will  give  the  current  density  J.  The  singularity  as  y  — ►  0 
is  best  treated  by  rewriting  Bz  as 
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where  t'  —  t  —  y2y.0(r/4s2.  The  form  of  (9)  is  well  behaved  as  y  —*  0.  The  current  is 
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Jy(x.y,t)  — - y  /  dk  ke~k  y  I\r)e~u  ,  (11) 

*  Jo 

where  u  —  k  [i-®0(t)]  ,  and  r  =  t  —  y,ocr/4k2.  The  simplest  special  case  is  where  the  current 
is  constant,  I'(t)  =  and  the  velocity  of  the  thin  armature  is  also  constant,  xo(t)  —  vt. 
The  above  integrals  for  the  current  can  then  be  evaluated  in  closed  form  throughout  the 
rails, 

Jx  =  ^e-^/2[tfo(fcor/2)  -  ^(feor/2)],  (12) 

Z7T  r 

and 

Jy  =  -e_fc°^^2/ir1(A:0r/2),  (13) 

Z7T  r 

where  fc0  =  /x0eru  is  an  inverse  length  that  can  be  used  to  scale  all  of  the  variables, 
r  -  Vt2  +  y2  is  the  distance  from  the  field  point  to  the  present  position  of  the  armature,  £ 
was  defined  below  (2),  and  the  K’s  are  Bessel  functions  of  the  second  kind  with  imaginary 
argument.  At  a  distance  from  the  current  source  large  compared  to  k^1,  the  current 
density  of  Eqs.  (12)  and  (13)  has  the  asymptotic  form8 


(14) 


From  here  one  can  show  that  the  shape  of  these  flow  lines  is  parabolic  at  large  distances 
from  the  origin.  In  Eq.  (14)  the  space  in  front  of  the  armature  has  £  >  0  and  r  >  0  so  the 
exponential  is  damped  very  quickly  (for  large  ko )  giving  negligible  current  ahead  of  the 
armature  as  expected.  Behind  the  armature,  £  <  o,  there  is  damping  for  y  >  0,  but  along 
the  rail  surface  the  exponent  is  zero. 

The  current  of  Eqs.  (10)  and  (11)  can  be  evaluated  asymptotically  for  a  general 
armature  current  and  velocity.  The  integrals  are  approximated  using  the  method  of 
steepest  descent;  We  have  done  this  calculation.  It  is  straightforward  but  somewhat  long. 
In  the  region  where  x0(t)  —  x  y  >  0  there  is  a  much  simpler  way  to  obtain  the  current. 
For  large  cr  the  argument  of  the  error  function  in  Eq.  (9)  varies  rapidly.  We  make  the 
approximation  that  the  error  function  is  a  step  function,  with  1  —  erf(x)  ~  2  for  x  <  0, 
and  1  —  erf(x)  m  0  for  x  >  0.  Then  the  magnetic  field  of  Eq.  (9),  behind  the  armature 
where  £  =  x  —  x0(t)<0,is  approximately 


Bx{x,y,t)  « f  dse  *7 l'(t  -  ay0y2/4s2), 
V*  J 3. 


(15) 
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with  sx  =  y/8 ,  and  we  have  introduced  a  skin  depth  8  defined  by 

«  =  (16) 

v  Ho<r 

where  iz  <  t  is  the  time  when  the  armature  passes  position  x ,  defined  by  x0(tx)  =  x. 
Eq.  (15)  has  a  pleasing  form  and  is  similar  to  some  of  the  results  in  Knoepfel.6  We  obtain 
the  currents  from  the  curl  of  Eq.  (15), 


J*(x,y,i)  ~  I'(t)e  y  /6 


Ju  di 
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Jy(x,y,t)  *  -=-  y  J'(tx)e~y  ,  (18) 

0 \J Ti  V x\t  *z) 

where  vx  is  the  velocity  at  time  tx.  We  have  written  Jx  in  Eq.  (17)  so  that  the  first  term 
contains  I'(t).  If  that  term  is  combined  with  the  second  term  in  the  integral  it  can  be 
seen  that  the  first  term  would  then  contains  V{tx).  We  note  that  even  if  I'(t)  remains 
positive,  Jx(x,y,t)  may  become  negative  because  of  the  term  with  dl'/dt ,  whenever  the 
current  decreases.  Eqs.  (17)  and  (18)  are  very  accurate  when  we  are  just  a  few  lengths 
(/ro<rv)_1  away  from  the  present  position  of  the  armature.  This  was  verified  by  extensive 
comparisons  with  numerical  evaluations  of  the  exact  expressions  in  Eqs.  (10)  and  (11). 

The  length  k^1  is  very  small  in  typical  rail  launcher  situations.  For  the  case  of  copper 
at  room  temperature,  and  a  speed  of  1  km/sec,  the  length  k^1  =  1.25  x  10~5  meters. 
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A  picture  of  the  current  flow  for  these  parameters  is  shown  in  Fig.  2.  The  acceleration  is 
8  x  107  m/s2.  Here  the  Hues  are  tangent  to  the  current  density  vector  J  and  the  curves 
are  drawn  at  equal  increments  of  current;  the  last  curve  includes  95%  of  the  total  current. 
The  current  has  passed  its  peak  and  is  decreasing.  Some  reverse  current  flow  is  evident  in 
the  graph  in  regions  where  Jx(x,y,t )  has  become  negative.  This  should  not  be  surprising. 
This  reverse  current  is  just  an  eddy  current  that  is  trying  to  keep  the  magnetic  field  in  the 
rail  from  dropping  in  value  as  the  total  rail  current  decreases.  The  total  horizontal  scale 
is  1.75  cm;  the  total  vertical  scale  is  expanded  to  0.12  cm. 

Breech  and  Muzzle  Voltages 

There  are  several  circuit  mode^s^'Used  f°  describe  rail  launchers  as  circuit  elements; 
these  require  expressions  for  the  voltages  in  terms  of  the  currents.  These  expressions  can 
be  obtained  by  application  of  the  theorem 

jt  J  B  •  dS  =  J  ^  -  V  x  (u  x  B)  +  (V  •  B)u  •  dS, 

where  the  surface  integral  is  taken  over  a  surface  moving  with  local  velocity  u.  For  the 
magnetic  field  B,  using  Maxwell’s  equations,  we  obtain 


J  B'dS  =  -  ^[E  +  ux  B]-  dr, 


(19) 


where  the  closed  line  integral  is  taken  over  the  moving  path.  In  a  moving  conductor  Ohm’s 
law  is  J  =  <t( E  -f  u  X  B). 

For  simplicity  we  will  derive  the  circuit  equations  using  the  thin  armature  results. 
Integrating  Eq.  (19)  around  the  path  3-4-5-6  shown  ahead  of  the  armature  in  Fig.  1,  we 
find  that  the  muzzle  voltage  Vm  is 


V5  -  v6  =  Vm(t)  =  Va(t)  =  J° 


Ja(t) 


dy34, 


(20) 


where  Ja  and  aa  are  the  current  and  conductivity  of  the  armature,  w  is  the  rail  separation, 
and  Va  is  the  resistive  voltage  drop  across  the  armature.  We  have  used  the  fact  that  B 
is  zero  ahead  of  the  armature  in  this  simple  geometry  and  that  the  current  is  extremely 
small  ahead  of  the  armature.  This  is  why  the  muzzle  voltage  in  Eq.  (20)  is  given  by  just 
the  resistive  drop  in  the  armature. 

We  use  the  path  of  integration  1-2-3-4  shown  behind  the  armature  in  Fig.  1  to  find 
that  the  breech  voltage  Vb  is 

Vb(t)  —  V2-V1 

d  r* o(t)  7  /y'l 

=  —  (B(t)wxo(t))  +  2  /  — — -dxu  +  Va(t),  (21) 

Clt  J  q  (T 

where  x0(t)  is  the  position  of  the  armature  at  time  t.  Here  we  have  used  the  fact  that  in 
this  geometry  Bz  is  uniform  behind  the  armature  and  that  the  integral  2  — >  3  is  the  same 
as  the  integral  4  — >  1. 


The  breech  voltage  Vb(t)  involves  the  integral  of  Jx  along  the  rail’s  inner  edge,  y  =  0. 
Equation  (17)  shows  that  Vb(t)  does  not  depend  on  only  the  values  of  /'  and  dl'  j dt  at  time 
t,  but  on  their  time  history.  The  circuit  equation  will  therefore  be  an  integro-differential 
equation.  If  dl1  /dt'  in  the  integral  is  not  too  large,  as  is  the  case  in  most  rail  launchers,  we 
can  simplify  the  result  by  doing  a  Taylor  expansion  around  the  upper  limit  of  integration 
and  carrying  out  the  t1  integral.  The  result  is 


J i(£,0,  t) 
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+  30^~tx)  dt 3 


which  is  an  expansion  in  powers  of  ( t  —  tx).  Especially  nearer  the  armature,  the  contri¬ 
butions  of  the  higher  derivatives  of  I'  in  Eq.  (22)  may  be  neglected  compared  to  the  first 
derivative.  We  can  now  use  these  results  in  writing  the  breech  voltage  of  Eq.  (21).  We  use 
the  fact  that  behind  the  armature,  between  the  rails, 

B(x,0,t)  =  fi0I(t)/ht 

where  h  is  the  rail  height,  and  I(t)  is  the  total  rail  current,  I{t)  =  hl'(t),  to  write 
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The  first  term  in  Vj, (f)  allows  us  in  this  simple  geometry  to  calculate  the  usual  rail 
inductance 

w 

L0{t)  w  (24) 

proportional  to  the  distance  x  that  the  armature  has  traveled.  The  second  term  is  the 
rail  resistance  term.  The  skin  depth  that  enters  there  involves  the  time  t  —  tx,  which  is 
the  time  since  the  armature  passed  position  x.  The  term  Va{t )  is  the  resistive  drop  in  the 
armature.  The  last  derivative  term  in  Eq.  (23)  is  a  skin  inductance  term  that  is  small 
compared  to  the  main  inductance  in  typical  railgun  situations.  The  ellipsis  indicate  the 
presence  of  second  and  higher  derivatives  from  Eq.  (22).  These  would  be  negligible  only 
for  slowly  varying  currents. 

The  effective  resistance  of  each  rail  can  be  read  from  Eqs.  (16)  and  (23).  We  use  the 
definitions  of  tx  and  8  to  write  the  rail  resistance  at  time  t  as 


m  =  - 


,  vo(t') 
Vt^V' 


where  vo(t)  =  dxo(t)/dt.  R(t)  can  be  expressed  in  terms  of  the  instantaneous  armature 
position  x  only  if  x  is  known  as  a  function  of  t.  We  note  that  an  elementary  derivation  of 
R(t)  would  miss  the  s/tk  factor.  For  example  if  the  armature  speed  is  constant  we  get 
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where  we  used  x  =  n0h  This  implies  that  there  is  no  fixed  resistance  per  length,  but  rather 
that  the  resistance  varies  more  slowly  with  length  than  linearly.  The  case  of  constant 
acceleration,  x  =  at2  j 2  gives 


**>  =  3WS“,/4(2*)I/4- 


Temperature  Distribution 


The  local  heating  of  the  rail  is  proportional  to  the  square  of  the  current  densi  ty.  Solve 
the  heat  diffusion  equation 

dT  ,  J2 
pc— -xv2r=-, 


neglecting  the  latent  heat  of  melting,  where  c  is  the  specific  heat  per  mass,  p  is  the  mass 
density,  K  is  the  heat  conductivity,  and  T  is  the  temperature.  The  general  solution  of  this 
equation  in  our  geometry  can  be  given  in  terms  of  Green’s  functions.  We  only  need  however, 
the  case  of  small  dimensionless  diffusion  coefficient,  D  =  Kp0a/pc  <C  1.  ( D  =  0.0053  for 
copper  at  room  temperature.)  In  this  case  one  can  verify  that  the  the  heat  diffusion  is 
negligible.  During  the  short  time  that  the  current  flows,  the  heat  generated  at  a  point 
in  the  rails,  /  dt  J2 / a ,  does  not  have  time  to  diffuse  away  but  rather  it  stays  where  it  is 
produced.  The  temperature  rise  at  a  point  is  then 


T(x,y,t)  =  —  f  dt1  J2(x,y,t') 
pac  J 


This  formula  implies  that  the  rate  of  heating  is  greater  at  points  near  the  arc,  but  it  also 
says  that  the  higher  temperatures  are  reached  at  points  farther  behind  the  arc  (closer  to 
the  breech).  This  happens  because  the  heat  diffusion  in  this  short  time  is  so  small,  and 
because  points  near  the  breech  are  subjected  to  the  current  for  a  longer  time.  Of  course 
radiation  cooling  can  affect  this. 

First  we  will  calculate  the  temperature  at  the  rail  surface  for  the  case  of  a  thick 
armature  of  length  2 L.  The  special  case  where  both  the  rail  current  and  the  velocity  are 
constant  in  time  allows  us  to  express  the  current  at  the  rail  surface  in  closed  form.  Due 
to  the  lack  of  heat  diffusion  this  will  be  the  most  important  part  of  the  current  needed  to 
understand  the  ohmic  heating  and  melting  of  the  surface.  The  y-component  of  the  surface 
current  density  is  zero  except  in  the  armature  region  from  £  =  —  L  to  £  =  +L,  where  it 
has  the  value  I' /2L.  Recall  that  |£|  is  the  distance  along  the  rail  from  the  point  with 
coordinate  x  to  the  center  of  the  armature.  The  i-component  of  the  current  density  at 
the  y  —  0  rail  surface  is  given  by 


J*(x,0,t)  =  ^  [M(M£  +  L))  -  M{hU  -  L)j 


where 


M(x)  -  {(x  +  l)K0(\x\/2)  -  |x|JsT1(|a:|/2)]e- 


The  temperature  at  the  surface  can  be  found  in  the  approximation  of  Eq.  (27)  using  the 
current  of  Eq.  (28).  The  integral  is  readily  evaluated  numerically.  The  temperature  is 
Ho I'~ /pc  times  a  dimensionless  factor  depending  on  k0L  and  &o£.  For  ko L  =  25  and 

=  —50000,  the  factor  is  3.04;  for  k0L  =  25  and  =  —5000  the  factor  is  2.31.  For 

k0L  =  2000  and  the  same  £’s,  the  factors  are  1.63  and  0.89  respectively. 

Again  we  take  copper  at  room  temperature  and  a  speed  of  1  km/sec,  where  the  scaling 
factor  ko  —  8  x  104  m-1,  so  a  ko(  of  50000  becomes  0.63  m;  koL  —  25  gives  an  armature 
width  2 L  of  0.6  mm;  koL  =  2000  corresponds  to  the  width  2 L  of  5  cm.  Using  the  room 
temperature  values  of  p  and  c,  and  a  typical  value  of  the  current  density  3x  107  Amp/ meter, 
the  factor  3.04  gives  a  temperature  rise  of  985  C.  The  factor  1.63  gives  528  C.  The  melting 
point  of  copper  is  1083  C  so  these  examples  do  not  lead  to  melting. 

The  temperature  rise  is  far  less  sensitive  to  the  length  of  the  arc  than  it  is  to  the  overall 
current  per  height  of  the  rail.  This  is  essentially  because  the  penetration  of  the  current 
into  the  rail  is  small  for  a  long  time;  so,  even  if  the  length  of  the  arc  is  long,  all  of  the 
current  will  eventually  have  to  pass  near  any  given  part  of  the  surface.  That  the  current 
does  gradually  move  into  the  rail  is  reflected  in  the  weak  sensitivity  of  the  temperature  to 
the  arc  length.  A  factor  of  80  in  koL  causes  a  factor  of  less  than  2  in  the  final  temperature. 
The  weak  sensitivity  to  the  armature  length  L  and  the  strong  I'2  dependence  implies  that 
the  arc  height  is  more  important  than  the  arc  length  in  raising  the  temperature.  Local 
pinching  of  the  arc  and  concentration  of  the  current  can  then  give  local  melting.  In  the 
examples  of  the  previous  paragraph,  the  copper  will  reach  its  melting  point  at  0.63  m  away 
from  a  6  mm  wide  armature  if  the  current  is  increased  by  5%.  In  the  case  of  the  5  cm  long 
armature,  a  40%  change  in  current  will  be  required. 

The  thin  armature  results  can  be  used  to  model  lateral  pinching  of  the  arc  into  a 
sheet  of  current.  In  reality  current  filaments  form  spots,  not  sheets  of  current.  However, 
many  filaments  moving  together,  as  are  sometimes  observed,  might  be  approximated  by  a 
sheet.  We  calculate  the  temperature  rise  in  the  case  of  the  thin  armature  with  constant 
current  and  velocity.  Remarkably  simple  analytic  results  will  be  obtained.  We  substitute 
the  asymptotic  form  of  the  current  in  Eq.  (15)  into  Eq.  (28)  and  obtain 
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where  r'  =  \/£'2  +  y2.  We  have  changed  variable  of  integration  in  Eq.  (27)  from  t'  to 
=  x  —  xo(t'),  with  xo(t )  =  vt,  and  we  have  let  I'  =  I',  where  the  subscript  s  indicates 
that  this  is  the  current  per  height  through  a  sheet.  The  upper  limit  of  zero  comes  from 
the  fact  that  there  is  very  little  current  ahead  of  the  armature,  so  a  point  begins  to  receive 
heat  essentially  only  after  the  armature  passes  it,  for  t'  >  tx,  or  £  <  0.  For  the  domain 
where  the  current  has  been  passing  for  a  time,  that  is,  for  £  <  0  and  |&0£[  1,  with  y  fixed 

and  small,  the  integral  can  be  approximated  still  further.  As  £'  varies  over  its  domain, 
the  quantity  +  r'  stays  nearly  equal  to  zero,  and  all  the  exponentials  are  then  equal  to 
one.  Deviation  from  this  approximation  occurs  only  in  the  neighborhood  of  £'  =  0.  The 
integrand  can  now  be  approximated  by  replacing  the  exponential  with  a  step  function. 


The  resulting  integral  is  then 


*PC  Ji  y/yl  +  t'2 

This  is  easily  evaluated  to  be 


irpc  \yj  ir  pc 


(30) 


From  Eq.  (30)  we  see  that  the  isotherm  with  temperature  Tm  is  the  straight  line  with 
equation 


For  copper,  the  coefficient  in  this  expression  is  about 
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The  remarkably  simple  expression  in  equation  (31)  has  several  interesting  conse¬ 
quences.  The  amount  of  surface  melting  can  be  estimated  from  it  if  we  neglect  the  latent 
heat  of  melting.  If  the  temperature  Tm  is  set  equal  to  the  melting  point  of  the  rail,  we 
can  find  the  depth  to  which  a  given  current  sheet  will  melt  the  surface.  For  copper,  and  a 
current  sheet  I'  =  300,000  Amps/cm  of  rail  height,  the  slope,  \y/x\,  of  the  melting  curve 
is  found  to  be  3.0  x  10“5.  The  mass  of  material  melted  will  then  be 


m  = 


-pR2hexp 


(32) 


where  p  is  the  density,  R  is  the  length  of  the  rail,  and  h  is  the  rail  height.  The  exponential 
dependence  of  the  melted  mass,  m,  on  the  rail  parameters  should  be  noted.  Taking 
/'  =  300,000  Amps/cm,  p  =  8.9  gm/cm3,  R  =  400  cm,  h  =  1  cm,  and  Tm  —  1,083  °C,  we 
obtain  m  —  5.2  grams  for  copper  rails.  The  large  amount  of  melting  is  due,  of  course,  to 
the  laxge  value  used  for  the  current  sheet. 

Conclusion 

The  calculations  above  have  been  done  for  an  idealized  geometry  where  there  are  no 
variations  in  the  z  direction.  These  simplifications  were  introduced  in  order  to  be  able  to 
treat  the  time  dependent  problem  analytically.  The  formulas  obtained,  however,  should 
be  of  value  in  understanding  the  performance  of  railguns.  We  believe  that  despite  the 
simplifications,  reasonable  estimates  of  rail  currents  and  temperatures  can  be  made  using 
our  methods. 

We  emphasize  one  of  the  most  interesting  results  of  this  paper.  This  is  the  possibility  of 
local  current  density  reversals  when  d,V fdt  becomes  negative.  This  is  simply  an  inductance 
effect.  The  current  reverses  direction  to  try  to  prevent  a  decrease  in  the  value  of  Bz(x,y,t ) 


in  the  rail.  The  same  reversal  can  occur  in  parts  of  the  armature  in  a  real  case.  The 
reversal  in  the  armature  does  not  occur  in  our  model  because  we  assume  a  known  current 
density  distribution  in  the  armature  and  rigorously  compute  the  current  in  the  rails.  In 
the  real  problem  the  current  distribution  in  the  ^rrnature  is  not  specified;  only  the  toted 
current  is.  Then  a  decrease  in  the  total  current  would  decrease  the  magnetic  field  in  the 
interior  of  the  armature.  The  current  density  in  the  rear  of  the  armature  could  easily 
reverse  then,  just  as  it  did  near  the  rail  edge  in  Fig.  2,  in  order  to  try  to  keep  up  the 
value  of  the  magnetic  field  in  the  interior  of  the  armature.  The  portion  of  the  armature 
where  the  current  is  reversed  could  then  be  subject  to  a  magnetic  force  directed  toward  the 
breech.  This  would  have  a  powerfully  disruptive  effect  in  the  case  of  a  plasma  armature. 
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Abstract 

We  report  on  our  development  of  a  two  dimensional  MHD  code  to  simulate  internal 
motions  in  a  railgun  plasma  armature.  We  use  the  equations  of  resistive  MHD,  with  Ohmic 
heating,  and  radiation  heat  transport.  We  use  a  Flux  Corrected  Transport  code  to  advance 
all  quantities  in  time.  Preliminary  runs  show  the  growth  and  subsequent  shedding  of  vortex 
structures  in  response  to  a  small  perturbation  upon  an  initial  equilibrium.* 

Introduction 

In  a  railgun  the  armature  carries  the  current  across  from  one  rail  to  the  other  and 
is  subjected  to  a  powerful  J  x  B  magnetic  force  that  pushes  it  against  the  rear  of  a 
projectile.  In  many  cases  the  strong  heat  produced  by  the  current  generates  a  plasma 
armature.  The  plasma  armature  is  then  a  hot,  high  density,  high  pressure  arc  undergoing 
rapid  acceleration.  The  magnetic  force  compresses  the  plasma  and  causes  a  high  pressure 
at  the  front,  where  the  plasma  pushes  against  the  projectile.  The  pressure  takes  a  profile 
that  diminishes  toward  the  rear  as  the  plasma  approaches  a  region  of  much  lower  pressure 
and  density  . 

Much  of  the  modeling  of  plasma  armatures  in  railguns  begins  with  the  one  and  two 
dimensional,  infinite  rail  height,  MHD  models  of  Powell  and  Batteh1,2.  Their  models 
describe  a  great  deal  of  the  physics,  including  radiation  transport,  and  the  degree  of  plasma 
ionization.  Those  models  assume  that  the  plasma  is  in  a  steady  state  with  no  flow  in  the 
frame  of  the  ac  derating  armature,  and  yield  profiles  of  the  magnetic  field,  and  of  the 
fluid  properties  in  the  armature.  Some  work  has  been  reported  by  Sloan3  on  extending 
the  steady  models  to  three  dimensions. 

The  steady  state  assumption  is  an  important  limitation  of  those  models.  The  calcu¬ 
lations  of  Huerta  and  Decker4'5,  as  well  as  those  of  Powell6,  give  that  the  steady  plasma 
armature  is  subject  to  the  MHD  flute  instability  that  occurs  when  a  magnetic  field  acceler¬ 
ates  a  plasma.  The  flute  instability  causes  two  dimensional  corrugations  of  the  rear  of  the 
plasma  that,  together  with  ensuing  flows,  may  destroy  the  simple  steady  states  found  by 
the  methods  of  Refs.  [1]— [3] .  Recently1’8,  studies  have  been  reported  on  one  dimensional 
time  dependent  models.  These  numerical  studies  emphasize  armature  initiation  phenom¬ 
ena,  and  are  not  able  to  assess  the  effects  of  two  dimensional  instabilities  on  the  overall 


steady  state. 

We  report  our  work  on  a  two  dimensional,  MHD,  infinitely  tall  rail,  time  dependent 
numerical  simulation  of  the  behavior  of  the  plasma  armature.  Our  model  is  capable  of 
describing  not  only  the  growth  of  the  flute  instability,  but  also  other  two  dimensional 


*  This  work  supported  in  part  by  the  Air  Force  Office  of  Scientific  Research  under  grant 
number  84-0116. 
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effects,  such  as  the  shedding  of  plasma  at  the  rear  of  the  arc,  and  convective  transport  in 
the  armature.  We  also  seek  to  describe  current  filamentation,  possible  current  reversal  due 
to  a  drop  in  total  current,  armature  disruption  due  to  the  rearward  forces  on  the  reversed 
currents,  as  well  as  a  variety  of  waves,  and  shocks. 

The  numerical  solution  of  the  governing  two  dimensional  MHD  equations  is  accom¬ 
plished  with  an  algorithm  which  is  designed  to  take  advantage  of  the  architecture  of  Multi- 
Vector  processor  supercomputers. The  basic  algorithm  for  advancing  the  critical  quantities 
in  time  is  a  finite  difference,  explicit,  Eulerian,  Flux  Corrected  Transport  (FCT)  sub-code 
which  closely  follows  an  algorithm  developed  at  NRL9,10. 

The  Governing  MHD  Equations 

The  plasma  is  collision  dominated  and  can  be  described  by  the  equations  of  magne¬ 
tohydrodynamics  subject  to  the  usual  limitations  of  MHD11.  The  plasma  is  treated  as  a 
single  fluid  of  density  p,  velocity  v,  pressure  p,  energy  per  unit  mass  e,  and  temperature 
T .  We  choose  the  x-axis  along  the  rail  direction,  from  the  projectile  back  toward  the 
breech,  and  the  y-axis  is  transverse  to  it,  opposite  to  the  direction  of  the  current  in  the 
armature,  with  the  z  axis  along  the  direction  of  the  magnetic  field  between  the  rails.  The 
projectile  moves  relative  to  the  rails  with  a  velocity  vp  toward  the  muzzle.  Our  coordinate 
system  is  fixed  to  the  projectile.  When  we  speak  of  the  fluid  velocity  v,  the  electric  field 
E,  or  any  other  quantity,  we  mean  it  as  measured  in  the  frame  of  reference  that  moves 
with  the  armature,  unless  otherwise  specified.  It  is  easy  to  show  that  the  only  relevant 
( vp  <C  c )  change  in  the  fluid  and  Maxwell  equations  brought  about  by  our  choice  of  moving 
coordinate  system,  is  the  appearance  of  a  ’’gravity  term  g”  in  the  momentum  equation, 
where 

dvP~  /i\ 

g  =  (1)° 

Therefore,  as  the  armature  gains  speed  toward  the  muzzle,  g  =  dvp/dt  >  0,and  g  points 
toward  the  breech.  All  quantities  are  taken  independent  of  z  in  our  two  dimensional 
model  (translational  invariance  along  the  z  axis,  0  =  d/dz  as  for  infinitely  tall  rails).  The 
magnetic  field  only  has  a  z  component,  but  the  fluid  velocity,  and  the  current  do  not  have 
z  components. 

In  MHD  the  relevant  Maxwell  equations  are 

V  ■  B  =  0,  VxE  =  -|,  (2)i 


V  x  B  =  p0J  .  (3)2 

The  displacement  current  has  been  neglected  in  Ampere’s  law,  Eq.  (3),  because  only 
frequencies  much  lower  than  cr/eo  are  considered.  The  electric  field  is  given  by  a  simple 
Ohm  s  law 

J  =  <r(E  +  v  x  B),  (4)? 

where  a  is  usually  taken  as  the  Spitzer  conductivity1.  Solving  Eq.  (4)  for  E,  eliminating 
J  using  Eq.  (3)  and  substituting  into  Eq.  (2)  one  obtains  the  equation  that  advances  B  in 


time 


dB 

dt 


(5)4 


-  +  V-(B,v)  =  —[■  1 


■-]  +  #:[- 


1  dB , 


]• 


dx  1  pocr  dx  J  dy  L/i0 cr  dy 

The  quantity  l//x0 cr  inEq.  (5)  is  a  coefficient  of  magnetic  diffusion.  When  the  magnetic 
diffusion  coefficient  is  small  (large  u),  the  magnetic  field  tends  not  to  diffuse,  but  rather 
to  be  convected  along  with  the  fluid.  When  the  diffusion  coefficient  is  large,  the  magnetic 
field  diffuses  through  the  conducting  fluid.  The  tall  rails  carry  a  current  I'  per  unit  height. 
A  real  rail  of  height  h  carrying  a  current  I  is  modelled  by  a  current 


h 


(6) 
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In  the  space  behind  the  armature,  where  there  is  no  current,  the  magnetic  field  is 
uniform  with  the  value 

Bz(t)  =  p0I'(t).  (7)4b 

In  MHD  the  plasma  is  treated  as  a  single  conducting  flilid  even  though  it  consists  of 
several  fluids.  The  inertia  is  carried  essentially  by  the  heavy  particles.  The  electron  fluid 
is  not  treated  explicitly  but  has  its  effect  via  the  current  J.  The  fluid  equations  are  the 
usual  ones.  First  the  equation  of  mass  conservation 


dp 
dt  + 


(?v)  =  0. 


(S> 


which  advances  p  in  time.  The  fluid  velocities  vx  and  vy  are  advanced  by  the  equations 
for  the  x  and  y  components  of  momentum, 


d(pvx)  ,  d  ,  1  dv „ 

+  V.(^,v)  —  B,)  +  p— , 


dt 


and 


^  +  V  ■  (ptvr)  ==  —  ±ip+±.Bl). 


(9)6 


(10)7 


In  Eqs.  (9)  and  (10),  the  magnetic  force  J  x  B  has  been  rewritten  by  eliminating  J  as 
usual  from  Eq.  (3).  This  makes  the  magnetic  force  appear  as  the  gradient  of  a  magnetic 
pressure  B2z/2po-  The  fluid  pressure  p  is  advanced  from  the  equation  for  the  energy  e  per 
unit  mass, 


— +  V  •  (pev)  +  pV  •  v  =  -V  •  q  +  —  +  — , 
Ctt  cr  a 


(11)8 


and  relating  p  to  e  using  an  energy  equation.  The  simplest  energy  equation  is  for  the 
monoatomic  ideal  gas  with  7  =  5/3, 


e  =  (7  -  1  )p/p  . 


(12). 


There  are  several  effects,  such  as  variations  in  the  degree  of  ionization, that  require  modifi¬ 
cation  of  Eq.  (12).  The  temperature  T  is  introduced  via  an  equation  of  state.  One  of  the 
simplest1  is  for  an  ideal  gas  of  ions  and  electrons 


p  =  (1  +  a)pkBT /mo, 


(13)io 


m 

.vSl 


*«» 
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where  a,  the  ratio  of  electrons  to  heavy  particles,  is  discussed  in  Ref.  [1],  and  m0  is  the 
mass  of  the  ion  or  neutral  particle.  The  current  dependent  terms  in  the  right  hand  side  of 
Eq.  (11)  represent  Ohmic  heating.  The  term  with  the  heat  flow  vector  q,  represents  heat 
transfer.  The  form  of  q  in  terms  of  T  for  radiation  transport  in  the  diffusion  approximation 
is  given  in  Ref.  [l] . 


Boundary  conditions 


The  front  boundary  of  the  plasma,  at  x  =  0  is  up  against  the  back  of  the  projectile. 
The  plasma  pressure  acting  over  the  rear  area  of  the  projectile  provides  the  force  that 
pushes  the  projectile  forward  and  must  satisfy  the  boundary  condition 


tw  dvP 

J  p(*  =  0,y,t)hdy  =  Tnp  —  , 


where  w  is  the  rail  separation. 

We  have  placed  a  low  density  region  of  zero  conductivity  inside  the  region  of  computa¬ 
tion  at  the  rear  of  the  plasma.  The  conductivity  is  set  to  zero  when  the  density  falls  below 
a  certain  threshold  value  which  we  take  approximately  five  hundred  times  less  than  the 
density  near  the  projectile.  This  boundary  between  a  region  of  conductivity  and  a  region 
of  no  conductivity  has  to  be  treated  carefully  since  it  is  not  a  fixed  boundary.  In  the  region 
of  no  current  the  fluid  properties  are  still  calculated  and  advanced  in  time  from  the  fluid 
equations.  The  magnetic  field,  however,  is  simply  uniform  in  the  rear  J  =  0  region,  and 
given  by  Eq.  (7).  As  soon  as  the  density  at  a  cell  rises  above  the  threshold  density  for 
conductivity,  its  properties  are  advanced  with  the  full  set  of  MHD  equations. 

The  rails  are  treated  as  infinitely  conducting.  We  take  the  boundary  condition  that  the 
tangential  electric  field  in  the  frame  of  the  rails  is  zero.  This  forces  the  normal  derivative  of 
the  magnetic  field,  dBz/dy  to  be  zero  at  the  rail  surfaces.  The  value  of  dBz/dy  at  the  rails 
is  needed  in  the  right  hand  side  of  Eq.  (5),  because  of  the  second  derivative  with  respect 
to  y,  to  advance  Bz  at  the  center  of  a  cell  adjacent  to  a  rail.  At  the  rear  of  the  projectile 
we  have  Jx  =  0,  which  sets  Bz  =  constant,  in  fact  zero,  along  the  x  —  0  boundary.  This 
enables  us  to  calculate  the  second  derivative  of  Bz  with  respect  to  x  at  a  cell  adjacent  to 
the  projectile.  The  other  boundary  conditions  are  standard  except  for  the  heat  flow  vector 
at  the  boundaries  of  the  plasma.  Here,  in  order  to  evaluate  V  •  q  in  the  right  hand  side 
of  Eq.  (11)  we  need  the  normal  component  of  q  at  each  boundary.  Physically  we  need 
a  statement  of  how  the  plasma  radiates  to  its  surroundings.  A  variety  of  conditions  are 
possible,  such  as  the  one  used  in  Ref.  [1], 


qn  =  2  crgT64 


where  qn  is  the  component  of  q  along  the  outward  normal  at  the  plasma  boundary,  and 
T),  is  the  plasma  temperature  at  the  boundary.  The  factor  of  2  in  Eq.  (15)  is  explained  by 
Zeldovich12.  The  radiation  that  leaves  the  plasma  actually  comes  from  some  distance  in 
its  interior,  where  the  temperature  is  a  bit  higher  than  at  the  surface,  21/4  higher. 


Numerical  Methods 


OCxTum  vx  vm  tfxuxvx 
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We  do  a  two  dimensional,  Eulerian  calculation  using  a  fixed  rectangular  array  of 
square  cells.  A  typical  grid  will  have  20  cells  in  the  y  direction  to  cover  the  one  cm  rail 
separation,  and  240  cells  in  the  x  direction.  The  main  conducting  plasma,  which  initially 
is  set  at  10  cm  length,  extends  from  cell  1  to  200.  Beyond  cell  200  we  have  a  low  density, 
nonconducting  gas  region  that  extends  to  the  end  of  the  computational  region  at  cell  240. 
As  the  flows  develop,  any  gas  on  a  cell  in  the  240th  row,  with  a  velocity  to  the  rear,  is 
allowed  to  escape. 

The  time  integration  uses  a  simple  midpoint  rule  which  is  second  order  accurate. 
The  driving  variables  are  determined  at  the  half  step  by  a  forward  differenced  first-order 
algorithm,  and  are  then  used  to  advance  all  the  variables  a  full  time  step  in  a  time  centered 
manner.  We  use  a  one  dimensional  FCT9,10  code  that  is  optimized  for  a  vector  processor. 
In  order  to  preserve  the  optimization  we  do  the  two  dimensional  problem  with  the  one 
dimensional  FCT,  and  a  time  step  splitting  technique13.  The  time  step  splitting  technique 
can  also  take  advantage  of  the  capabilities  of  multi-processor  machines. 

This  FCT  is  an  explicit,  finite  difference  technique  that  is  essentially  second  order 
accurate,  and  has  strict  limits  on  its  stability.  It  imposes,  however,  the  further  conditions 
that  densities,  such  as  p,  and  e  remain  positive.  This  makes  it  of  indeterminate  order,  but 
enables  it  to  yield  physically  realistic  and  accurate  results,  even  in  the  presence  of  inviscid 
shocks  and  steep  gradients.  The  algorithm  provides  fourth  order  accurate  phases  and 
minimum  residual  diffusion.  It  consists  conceptually  of  two  stages.  The  usual  transport 
stage  is  followed  by  a  corrective  antidiffusive  stage.  FCT  has  been  used  by  Emery14,  for 
example,  in  the  study  of  Rayleigh-Taylor  instabilities  on  the  surface  of  targets  accelerated 
by  laser  ablation. 


Results 

The  code  is  optimized  to  run  on  a  supercomputer  with  vector  processing,  but  so  far 
we  have  only  run  preliminary  tests  on  a  VAX  8650,  where  it  runs  rather  slowly.  With  a 
grid  of  20  x  240  cells,  it  advances  some  1,000  time  steps  in  about  one  hour  of  CPU  time. 

We  report  preliminary  results.  We  have  set  up  programs  to  do  simulations  with 
parameters  similar  to  Ref.  [1].  A  typical  case  has  a  plasma  of  copper  ions,  mo  =  1.1  x  10~25 
kg,  with  a  projectile  mass  of  mp  =  2  gm,  and  a  plasma  mass  of  .065  gm.  The  railgun  cross 
section  has  h  =  1  cm,  and  w  =  1  cm.  The  preliminary  runs  have  an  adiabatic  energy 
equation  (zero  right  hand  side  in  Eq.  (11)),  and  all  the  atoms  doubly  ionized,  with  degrees 
of  ionization  xi  =  0,X2  =  1,  as  defined  in  Ref.  [1].  For  the  current  per  unit  height  we  take 
I'  =  2  x  107  Amp/m.  This  produces  a  magnetic  field  of  8tt  T  from  Eq.  (7).  The  pressure 
on  the  projectile  is  of  the  order  of  2,500  atmospheres  and  the  projectile  acceleration  is 
of  the  order  of  1.2  x  IQ7  m/sec2.  The  preliminary  runs  have  a  constant  conductivity 
a  =  1  x  105  (ohm-m)_1,  that  is  turned  off  to  zero  where  the  plasma  density  falls  below 
2.61  x  10~2kg/m2  at  the  plasma  rear.  The  linear  growth  rate  of  the  flute,  or  interchange 
instability,  is  about 

7  =  yfkg  =  (16)i* 

Eq.  (16)  is  the  same4,6  as  the  Rayleigh-Taylor,  or  Kruskal-Schwarzschild  rate.  According 
to  this,  we  expect  the  interchange  instability  to  develop  with  a  growth  rate  of  8.6  X  104 


The  grid  has  240  cells  along  the  x  axis.  At  t  =  0  we  set  up  a  plasma  equilibrium 
that,  except  for  the  temperature,  is  pretty  much  as  calculated  in  Ref.  [1].  The  plasma 
is  10  cm  long  and  extends  to  cell  200.  The  initial  temperature  is  made  uniform  and  set 
at  6.72  x  104  K.  The  plasma  is  initially  at  rest  in  the  frame  of  the  projectile  except  for  a 
velocity  perturbation  of  10  m/s  in  vx  from  cell  199  to  201.  The  velocity  perturbation  varies 
sinusoidally  in  y  with  a  wavelength  of  2/3  cm  or  about  13  y-cells.  The  vx  perturbation 
is  set  positive  (rearward)  near  the  walls,  and  negative  near  the  center.  In  the  same  cells 
we  also  perturb  the  pressure  by  2%  of  its  local  equilibrium  value  with  the  same  variation 
in  y  as  the  velocity  perturbation.  We  also  perturb  the  pressure,  again  by  2%  of  its  local 
equilibrium  value,  in  the  interior  of  the  plasma,  corresponding  to  cells  99-101. 


We  find  that  the  calculation  is  numerically  unstable  for  time  steps  longer  than  about 
At  —  2  x  10_9sec.  This  is  as  expected  from  the  Courant15  condition,  <  Ax,  due  to 

the  large  local  velocities  produced  in  the  low  density  region.  To  simulate  a  shot  down  a 
one  meter  long  barrel  requires  a  flight  time  of  about  0.4  X  10-3  sec.  This  amounts  to  some 
2  x  105  time  steps.  In  the  VAX  8650  this  would  amount  to  some  200  hours  of  CPU  time. 


We  show  the  results  of  a  run  that  used  At  =  1.86  x  10~9  sec,  and  went  for  40,500 
steps,  a  total  time  T  =  7.5  x  10~5  sec,  and  required  about  40  hours  of  CPU  time  on  the 
VAX  8650.  According  to  the  linear  growth  rate  of  Eq.  (16),  the  perturbation  should  grow 
by  e'yT  %  e6-4.  Fig.  1  shows  the  pressure  profile  at  t  =  0.  The  initial  2%  perturbation  is 
practically  invisible  in  the  profile. 


Fig.  1.  Pressure  profile  at  t=0,  s=0,  cells  180-210. 


Fig.  2  shows  the  current  J  after  7,000  time  steps,  t  —  13  x  10~6sec,  distance  travelled 
s  =  1.0  x  10_3m,  projectile  velocity  vp  =  1.6  x  102m/sec.  and  projectile  acceleration 
op  =  1.22  x  107m/sec2.  It  shows  that  the  conducting  boundary  has  oscillated,  and  a 
dominant  mode  is  developing. 
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?ig.  2.  J  at  7,000  steps,  t  —  13  x  10  6  sec,  s  =  1.0  x  10  3m,  cells  180-210. 


Fig.  3  shows  the  pressure  profile  after  17,500  time  steps,  t  =  3.3  X  10  5sec,  s  - 
6.5  x  10-3m,  vp  =  4.0  x  102m/sec  and  ap  =  1.22  x  107m/sec2.  It  has  a  noticeable  bulge. 


Fig.  3.  Pressure  profile  at  t  =  3.3  x  10  5sec,  s  =  6.5  x  10  3m,  cells  180-210. 


Fig.  4  shows  the  pressure  profile  after  27,000  time  steps,  t  =  5.0  x  10  5sec,  s 
1.5  x  10_2m,  vp  —  6.1  x  102m/sec  and  ap  =  1.22  x  107m/sec2. 


Fig.  5  shows  the  current  also  after  27,000  time  steps.  It  shows  further  boundary 
deformation. 


•  i  i 
i  i  ( 
.  i  i 
i  <  i 


I  I  t 
i  i  i 
i  i  : 


J 

i 


i 


•'  1  '  i 

1  1  1  i 

'  •  ‘  i 

1  1  1  i 

i  <  i  i 

tit'- 

i  i  / 

i  i  / 


1  i 
1  / 
i  / 


• .  '•ISF 


Fig.  5.  J  at  27,000  steps,  t  =  5  x  10  5  sec,  s  —  1.5  x  10  2m,  cells  180-210. 


Fig.  6  shows  the  velocity  field  v  also  after  27,000  time  steps,  from  cells  160-190,  not 
quite  to  the  end  of  the  conducting  region.  It  corresponds  to  the  pressure  profile  of  Fig.  4. 
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Fig.  6.  v  at  t  =  5.0  x  10  5sec,  5  =  1.5  x  10  2m,  cells  160-190. 


Fig.  7  shows  the  pressure  profile  after  40,500  time  steps,  t  =  7.5  x  10  5sec,  s 
3.5  x  10“2m,  vp  =  9.2  x  102m/sec  and  ap  =  1.22  x  107m/sec2. 


Fig.  7.  Pressure  profile  at  t  —  7.5  x  10  5sec,  s  =  3.5  x  10  2m,  cells  185-215. 


Fig.  8  shows  the  current  also  after  40,500  time  steps.  It  shows  that  the  boundary 
deformation  is  beginning  to  get  important.  Even  though  the  armature  has  travelled  only 
3.5  cm,  there  is  already  a  bit  of  plasma  shed  at  the  rear.  The  linear  growth  rate  gets 
larger  for  the  shorter  wavelength  modes.  However,  they  do  not  appear  to  grow  to  large 
amplitudes.  This  is  due  to  nonlinear  saturation  effects,  as  also  found  by  Emery14. 
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Fig.  8.  J  at  t  =  7.5  x  10  ssec,  s  =  3.5  x  10  2m,  cells  185-215. 
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